	PRO ICREC,N,fr,apert,kx1,ky1
; Initial conditions for magnetic reconnection 2.5 D
; Double precision numbers
; Input: N = # points, fr = width of the current sheet
; apert = amplitute of perturbation with kx1 and ky1 

B0 = 1
y = 2d*!dpi*dindgen(N)/(N)  ; periodic boundary conditions
x = 2d*!dpi*dindgen(N)/(N)  ; periodic boundary conditions
loadct,3 ; paleta de colores

bxdf = dblarr(N);
bx = dblarr(N,N) ;
a = dblarr(N,N);

C1 = -2d*fr*alog(cosh(.25d/fr)/cosh(.75d/fr))  ; make sure mean(a)=mean(bx)=0

for j = 0,N-1 do begin
	logo = alog(cosh((y(j)-1.5d*!dpi)/(2d*!dpi*fr))*cosh(.25d/fr)/(cosh((y(j)-.5d*!dpi)/(2d*!dpi*fr))*(cosh(.75d/fr))))
for i = 0,N-1 do begin
	pert = apert*cos(ky1*y(j))*cos(kx1*x(i))
	bx(i,j) = B0*(tanh((y(j)-1.5d*!dpi)/(2d*!dpi*fr))-tanh((y(j)-.5d*!dpi)/(2d*!dpi*fr))+C1)
	a(i,j) = B0*(2d*!dpi*fr*logo + C1*y(j)) + pert
endfor
endfor

bxdf(1:N-2)=(a(0,2:N-1)-a(0,0:N-3))/(y(2:N-1)-y(0:N-3)) ; bx from fin. diff.
bxdf(0)=(a(0,1)-a(0,N-1))/(y(2)-y(0)) ; first one 
bxdf(N-1)=(a(0,0)-a(0,N-2))/(y(2)-y(0)); last one 

window,2,xsize=300,ysize=300
plot,y,bxdf,linestyle=1,xr=[0,2.*!dpi],yr=[-2.,2.],xsty=1,ysty=1,xtit='y',ytit='Bx(x,y)',tit='Initital Conditions of B(x,y) using finite differences',col=black,back=1
oplot,y,bx[0,*]

window,1,xsize=300,ysize=300
plot,y,a[0,*],xr=[0,2.*!dpi],yr=[-2.,2.],xsty=1,ysty=1,xtit='y',ytit='a(x,y)',tit='Initital Conditions of a(x,y)',col=black,back=1

Nlev=10
alevp= 0 + (1+findgen(Nlev))*(max(a) - 0)/(1 + Nlev)
alevn=min(a)+(1+findgen(Nlev))*(0 - min(a))/(1 + Nlev)
window,3,xsize=N,ysize=N,title='a(x,y)'
contour,a,x,y,lev=alevp,xr=[0,2.*!dpi],yr=[0.,2.*!dpi],xsty=1,ysty=1,xtit='x',ytit='y',tit='Initial Countour Labels of a(x,y)',col=black,back=1
contour,a,x,y,lev=alevn,c_linest=1,/overplot
contour,a,x,y,lev=0,th=2,/overplot

;tvscl,a
;xr=[0,2.*!dpi],yr=[0.,2.*!dpi],xsty=1,ysty=1,xtit='x',ytit='y',tit='Initial Countour Labels of a(x,y)',col=black,back=1

; Calculo y grafico: J(x,y)
jj = dblarr(N,N)
k2 = dist(N)^2
ak = fft(a,-1)
jj = float(fft(k2*ak,1))
window,4,xsize=N,ysize=N,tit='j(x,y)'
tvscl,jj
;contour,jj,xr=[0,256],yr=[0,256],xsty=1,ysty=1,xtit='x',ytit='y',tit='Initial Countour Labels of j(x,y)',col=black,back=1

meanb = mean(bx)
meana = mean(a)
print,'Valor Medio a(x,y),B(x,y): ',meana,meanb
print,'a(x=0,y=0),a(x=0,y=2*pi) ',a(0,0),a(0,N-1)
print,'b(x=0,y=0),b(x=0,y=2*pi) ',bx(0,0),bx(0,N-1)
print,'bxdf(y=0),bxdf(y=2*pi) ',bxdf(0),bxdf(N-1)
end
